Polydisperse hard spheres at a hard wall 
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The structural properties of polydisperse hard spheres in the presence of a hard wall are inves- 
tigated via Monte Carlo simulation and density functional theory (DFT). Attention is focussed on 
the local density distribution p{a,z), measuring the number density of particles of diameter a at a. 
distance z from the wall. The form of p{a, z) is obtained for bulk volume fractions 77b = 0.2 and 
r]b = 0.4 for two choices of the bulk parent distribution: a top-hat form, which we study for degrees 
of polydispersity 5 — 11.5% and S = 40.4%, and a truncated Schulz form having S = 40.7%. Excel- 
lent overall agreement is found between the DFT and simulation results, particularly at rjb = 0.2. 
A detailed analysis of p{a, z) confirms the presence of oscillatory size segregation effects observed 
in a previous DFT study (Pagonabarraga et al., Phys. Rev. Lett. 84, 911 (2000)). For large 5, 
the character of these oscillation is observed to depend strongly on the shape of the parent distribu- 
tion. In the vicinity of the wall, attractive cr-dependent depletion interactions are found to greatly 
enhance the density of the largest particles. The local degree of polydispersity 5{z) is suppressed in 
this region, while further from the wall it exhibits oscillations. 



I. INTRODUCTION AND BACKGROUND 

Many complex fluids, whether natural or synthetic in 
origin, are intrinsically polydisperse in character, that 
is they comprise mixtures of similar rather than strictly 
identical constituent particles. Common examples are 
colloidal dispersions, macromolecules in solution and liq- 
uid crystals, where variation amongst the particles can 
occur in attributes such as their size, shape or surface 
charge. Polydispersity is of considerable practical rele- 
vance because it can affect the properties of materials in 
a variety of applications such as coating technologies [1] , 
foodstuffs [2] and polymer processing [3]. However, from 
the fundamental perspective, the current understanding 
of polydisperse fluids is considerably less advanced than 
that of their monodisperse counterparts. The reason for 
this is the inherent complexity that polydispersity im- 
parts to a system, and which generates difficulties for ex- 
perimental, theoretical and simulation approaches alike. 

On the experimental front, intricate technical issues 
beset the preparation, characterization and analysis of 
polydisperse samples, and only relatively recently has 
work begun to systematically address the generic conse- 
quences of polydispersity, such as its effects on phase be- 
haviour and the fractionation of particles of different sizes 
between coexisting phases [4-6] . As far as analytical the- 
ory is concerned, the principal challenge is the multitude 
of variables required to accurately describe the system's 
properties. Statistical mechanical theories of polydisper- 
sity typically regard the polydisperse attribute as a con- 
tinuous variable (tr, say) [7, 8]. Accordingly the system 
may be regarded as a mixture of an infinite number of 
species-each labelled by the value of cr. It is then natural 
to describe the polydispersity in terms of a density dis- 
tribution, p^{(t), measuring the number density of each 
species. Difficulties arise, however, when one attempts to 
determine the thermodynamic properties of the system, 
such as its phase behaviour. Because the free energy de- 



pends on (is a functional of) the entire distribution p^{cf), 
it occupies a parameter space that is effectively infinite 
dimensional. This in turn renders the analysis of phase 
behaviour much more problematic than for monodisperse 
systems or finite mixtures of a few components. 

For small degrees of polydispersity, perturbative ap- 
proaches permit some theoretical headway to be made 
[9, 10]. A more general approach is the moment free 
energy (MFE) method [11-14]. For certain truncatable 
forms of the free energy use of this method allows the 
effective dimensionality of the problem to be reduced to 
a manageable level by projecting the free energy onto a 
suitably defined subspace spanned by a few principal mo- 
ments of p°(cr). The phase behaviour of the resulting pro- 
jected free energy can be calculated relatively straight- 
forwardly. This method has been applied to investigate a 
range of open questions related to bulk phase equilibria, 
including (inter- alia) issues concerned with liquid- vapor 
coexistence and critical point shifts [15], freezing [16, 17] 
and liquid crystal phase transitions [18]. 

These advances in theoretical methods for dealing with 
polydisperse fiuids, have recently been complemented by 
parallel developments in simulation methodologies. New 
Monte Carlo (MC) algorithms now permit the effective 
study of polydisperse phase equilibria within the grand 
canonical ensemble (GCE). Use of this ensemble is ad- 
vantageous for the study of phase transitions because it 
allows the system density to fiuctuate as a whole, thereby 
catering naturally for both order parameter fluctuations 
and fractionation effects [21]. Within the GCE frame- 
work, polydispersity is controlled by means of a chem- 
ical potential distribution /i(cr), the form of which is 
adapted (for each state point of interest) in such a way 
as to yield an ensemble averaged density distribution 
p{a) that matches some prescribed form p^{cr)- This 
approach extends previous ones (see eg. refs.[19, 20]) 
by facilitating the targeting of a specific density distri- 
bution. It thus corresponds more closely to the experi- 
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mental situation for polydispcrsc fluids such as colloids 
and polymers, where the form of p{a) is fixed by the 
synthesis of the fluid and only its scale can vary depend- 
ing on the quantity of solvent present. The new tech- 
niques have recently been applied to obtain the equation 
of state of polydisperse hard spheres [21] and to inves- 
tigate polydispersity-induced alterations to bulk phase 
behaviour in a Lennard- Jones fluid [15]. 

While there is considerable ongoing progress in the 
study of bulk phase behaviour, to date comparatively 
little work has been reported regarding the influence of 
polydispersity on the properties of inhomogencous flu- 
ids [22]. One of the standard contemporary theoretical 
approaches for dealing with such systems is density func- 
tional theory (DFT) [23] . Recently, it has been shown by 
one of us, that the MFE method outlined above carries 
over straightforwardly to certain well established density 
functionals [24, 25] . The first such study of an inhomoge- 
ncous system considered the prototype model for a inho- 
mogencous polydisperse fluid, namely hard spheres in the 
vicinity of a hard wall [24]. Novel oscillatory behaviour 
was reported in the local concentration of different sized 
species as a function of the wall distance. Other DFT 
work has investigated polydisperse polymer /solvent mix- 
tures at coexistence [25], and the efi'ect of polydispersity 
on the fractionation and surface tension of a liquid-gas 
interface within a van dcr Waals approximation [26] . 

The extension of DFT techniques to polydisperse flu- 
ids is a welcome development. However, as yet, the relia- 
bility of the approximations inherent in such treatments 
remain untested. Clearly, therefore, it is desirable to ob- 
tain benchmarks with which these and other theories of 
confined polydisperse fluids can be compared. In this pa- 
per we address this issiie using MC simulation and DFT 
applied to size-disperse hard spheres at a hard wall. For 
this purpose the GCE simulation methods developed for 
the study of the bulk phase equilibria of polydisperse flu- 
ids carry over directly, allowing us to probe confinement- 
induced changes to the density distribution. In addition 
to facilitating a detailed comparison between simulation 
and DFT, our results confirm the presence of the oscilla- 
tory size segregation effcicts seen previously in the original 
DFT study, as well as further elucidating the character 
of the local fluid structure and the effects of depletion 
interactions near the wall. 

Our article is organized as follows: in Sec. II we outline 
the MC simulation scheme and the density functional 
theory. In Sec. Ill we analyse and compare in detail the 
MC and DFT results for two choices of the parent size 
distribution at two volume fractions. Finally, a discussion 
of the results and the prospects for interesting further 
work features in Sec. IV. 



II. SYSTEM AND METHODOLOGY 

The system with which we shall be concerned is a fluid 
of hard spheres, whose polydisperse attribute a corre- 



sponds to the sphere diameter. Pairs of particles i and j 
are assumed to interact via the potential 

TT t ^ / +«2 if |ri-r, | < 

C/,,(r„r,,a„a,) = | ^ ^^^^^.jj 2 . (i) 

Consider initially such a system in the bulk. Conven- 
tionally the state of the system is described by a "parent" 
density distribution [14, 27] P°{(t), which can be written 

= no/(a) . (2) 

Here no = N/V is the overall particle number density, 
while V is the system volume and /(cr) is a normalized 
shape function, the relative width of which is quantified 
by the dimensionless degree of polydispersity: 



measuring the standard deviation of the parent distribu- 
tion, normalized by its mean. 

Let us now break the translational symmetry of the 
system in one coordinate direction (z, say) by introduc- 
ing a hard wall in the plane z = 0. Interactions of the 
particles with the wall are assumed to be controlled by 
the potential: 

C/p.(^i,a,) = | otherwise " (4) 

Near such a wall, modifications to bulk behaviour arise 
even in monodisperse systems from standard packing ef- 
fects. For polydisperse systems, however, additional fac- 
tors come into play. Firstly, the excluded volume con- 
straint associated with Vpw implies that the centers of 
small particles are permitted to approach the wall more 
closely than those of large particles. Hence one expects 
(for sufficiently small z) that the local density distribu- 
tion will be truncated at large u with respect to the bulk 
(parent) form. Secondly, experience with binary hard 
sphere mixtures [28] suggest that very close to the wall 
one can expect attractive depletion interactions between 
the particles and the wall. In the polydisperse case these 
depletion forces should be cr-dependent. In order to quan- 
tify the net influence of all these effects, we consider the 
ensemble averaged local density distribution p(cr, z) at a 
perpendicular distance z from the wall: 

-2) = / / Pip^ r) dxdy . (5) 

In the present work, we have elected to study two forms 
of the parent distribution: a top-hat distribution and a 
Schulz, the form of which are shown in fig. 1. The top-hat 
distribution is defined by 
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fth{cr) 



(2c)-i ifl-c<(j<l + c 
otherwise 



(6) 



where (without loss of generahty), the mean particle di- 
ameter has been set to a = 1.0. We have studied the 
cases c = 0.2 and c = 0.7, for which 5 = c/\/3 « 0.115 
and 5 « 0.404 respectively. 

The Schulz distribution is defined by the shape func- 
tion 



1 (z + l 



z+l 



a exp 



(7) 



Here the parameter z controls the width of the distribu- 
tion and thence the value of S. We have considered the 
case z = 5, corresponding to S = {z + 1)~^/^ ~ 0.408. 
Again, by construction, a = 1.0. Notice that in con- 
trast to eg. a Gaussian, the Schulz distribution vanishes 
smoothly (has a natural cutoff) as cr ^ 0. For the pur- 
poses of the MC simulations described below, however, 
one does require an upper cutoff in a. Wc chose to trun- 
cate /(cr) at cTc = 4.0, which reduces the degree of poly- 
dispersity to 5 « 0.407. The truncated form of f{a) was 
then rescaled to unit integrated area. Given this trunca- 
tion, the population of particles at the cutoff diameter is 
very small compared to those having the mean diameter: 
p°(crc)//9°(a-) = 1.7 X 10~^. However, the volume of the 
largest allowed particle is 64 times that of one having the 
mean diameter. 
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FIG. 1: The top-hat and Schulz distributions considered in 
this work. 

A quantity useful in characterizing a polydisperse sys- 
tem is its volume fraction, which provides a measure of 
the degree to which the particles pack the available space. 
We shall employ it in preference to the number density 
no which does not alone provide equivalent information. 
We define the volume fraction in terms of the bulk parent 
distribution: 



% = J y da p°{a)a^. 



(8) 



For a nominated shape function /(cr), choosing a value for 
rjh (or equivalently no) serves to fix the bulk chemical po- 
tential distribution /i(cr). The assumption that the con- 
fined fluid exists in equilibrium with a bulk reservoir then 
implies that both the bulk and confined systems have 
equal ^(a). Note however, that for the confined hard 
sphere systems we consider here, the system-averaged 
volume fraction is generally less than that of the reser- 
voir. For both the top-hat and Schulz parents described 
above, we consider the cases % = 0.2 and % = 0.4. 

With regard to our use of scales, we shall throughout 
express all lengths in units of the average particle diam- 
eter a of the parent distribution. Owing to the lack of a 
temperature scale in hard sphere systems, we implicitly 
adopt the convention of assigning /? = l/ksT = 1. 



A. Grand-canonical ensemble Monte Carlo 
calculations 

The Grand Canonical Metropolis Monte Carlo algo- 
rithm we have employed has previously been described 
in ref. [21] and involves particle insertions, deletions, 
translations, and resizing operations. Particle resizing 
operations (a a'), or the insertion of a new particle 
are performed by drawing a random number with uni- 
form probability in the range of allowed sizes. Accep- 
tance/rejection criteria for these moves involve not only 
the change in the particle-particle and particle-wall po- 
tentials, but also the change in the chemical potentials 
IJ,{a) and n{a'). Translational moves (u, r) — » (cr, r') are 
not strictly necessary in the GCE formalism, as they are 
equivalent to the successive deletion of particle (cr, r) and 
insertion of particle (cr, r'), but their use proved benefi- 
cial to the sampling efficiency at the higher local volume 
fractions encountered near the wall. 

Within the GCE framework, the form of the ensem- 
ble averaged density distribution p{a) is controlled by its 
conjugate chemical potential distribution /x(cr). In order 
to perform simulations of a prescribed parent distribution 
one therefore needs to match p{a) to p^ia) in a bulk (or, 
in our case, fully periodic) simulation. Unfortunately, 
the task of determining the requisite /i(cr) is complicated 
by the fact that it is a functional of p^ijj). Effectively 
therefore, one is faced with solving the inverse problem 
/i(cr) = /i[p'^(cr)]. Doing so is facilitated by a recently 
proposed MC scheme-the non-equilibrium potential re- 
finement (NEPR) algorithm [29]- use of which enables 
the efficient iterative determination of ^[p°{a)], from a 
single simulation run, and without the need for an initial 
guess of its form. To achieve this, the method continually 
updates n{a) in such a way as to correct for the deviation 
of the instantaneous density distribution p{(j) from the 
target form (i.e. the parent). However, tuning /i(cr) in 
this manner clearly violates detailed balance. To counter 
this, successive iterations reduce the degree of modifica- 
tion applied to /i(cr), thereby driving the system towards 
equilibrium and ultimately yielding the equilibrium form 
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In order to obtain the bulk form of 12(a) correspond- 
ing to the top-hat and Schulz distributions at r]b = 0.2 
and rjb = 0.4, the NEPR algorithm was deployed in sim- 
ulations of a fully periodic cubic simulation cell of side 
L = 12. The resulting forms of fJ.{(j) for the respective 
parents were then used to study the effects of introduc- 
ing two oppositely facing hard walls at z = and z = L 
(the system remaining periodic in the x and y directions). 
The cell size L = 12 proved sufficient to ensure that the 
fluid properties at the cell midpoint were indistinguish- 
able (to within statistical uncertainty) from those of the 
fully periodic system. This finding was further tested by 
comparing the results with those for a cell that was elon- 
gated along the z direction, having (L^ = 24). Again 
within statistical uncertainties, no differences were ob- 
served, leading us to conclude that the two walls do not 
interact with one another and that the properties near 
one wall are representative of the semi-infinite geometry. 

The particle density p{(T,z) was accumulated in the 
form of a histogram in the course of the simulations and 
will be the object of our analysis in the next section. 
The histogram was formed by discrctising the permitted 
ranges < a < cTc and < z < L^, into bins. The 
bin width in a used for the top-hat distributions was 
5a = 0.01, while Sa = 0.04 was employed for the Schulz. 
The bin width in the z direction was Sz = 0.005 for 
the top-hat distributions and Sz = 0.02 for the Schulz 
distribution. 



B. Density Functional Theory model 

The density functional model employed is that intro- 
duced by Rosenfcld [23, 30], following the implementa- 
tion proposed in ref. [31], and corresponds to a weighted 
functional that has proved accurate for monodisperse sys- 
tems [32]. Although improvements have been proposed 
to the original Rosenfeld model which appear to perform 
better in the context of describing crystalline phases, we 
do not expect our implementation to suffer major draw- 
backs for the regimes of volume fraction regimes of inter- 
est in the present work. 

The free energy functional can be expressed as: 



n = y drda{p{a, r) [in {A^{a)p{a, r)) - l] 

+ [V{a,r) - ij{aMa,r)} + J dr.F^-[m„(r)] (9) 

where the excess contribution to the free energy func- 
tional, = — TOoln(l — TO3)-|-mim2/(l — m3)-|-m|/(l — 
ms)^, is a function of four moments only: 



ma(r) = / dadr' p{a,r)u)a{o;\r — r'\). (10) 



The four weight functions, LOa (non-local in space) are 
chosen to recover the Percus-Yevick free energy and cor- 
relation functions for a homogeneous mixture. In the 
bulk the moments have a simple physical interpretation. 
Indeed, 77x3(00) = rjb, mo(oo) = pb, while m2{oo)/pb is 
the mean area and mi(oo)/pf, the mean radius. 

The equilibrium density profiles are obtained by mini- 
mizing the free energy functional, eq.(9), with respect to 
the density profiles, leading to 

Pe,{a, r) = pO(<7) exp [-pV{a, r) + /.-(a, r) - n^a)] 

(11) 

where V{a,r) accounts for an external potential, and 

l^l'^ia) and p'^{a) correspond to the excess chemical po- 
tential and density of species a in the reservoir, respec- 
tively. 

The advantage of the moment structure of the func- 
tional is that the excess chemical potential of each species 
is 

drMa,\r-r'\)^ir') 
dvu:p{a,\v-v'\)p,f{v'), (12) 

13 ■' 

and can be interpreted as a linear combination of the four 
moment excess chemical potentials, = dT'^^ / drUa- 
This moment structure simplifies the study of polydis- 
persc mixtures [14]. 

For the case of a fluid mixture in the presence of a pla- 
nar hard wall, the external potential reduces to eq.(4): 
V{a,v) = Upw{z,a). The equilibrium density profiles, 
as given implicitly in eq. (11), are found numerically via 
iteration. To this end, it is necessary to deflne an under- 
lying lattice and to represent the parent distribution as a 
set of discrete species. However, the existence of this lat- 
tice imposes a minimum resolvable length scale (lattice 
cutoff). To avoid lattice artifacts, one must therefore ap- 
ply a lower cutoff to the parent distribution p^{a) which 
is itself large compared to the lattice spacing; we take the 
lower a cutoff to vary between 20 and 50 lattice units. 
Although in principle, this cutoff may lead to discrepan- 
cies between the DFT and MC results (in the simulations 
particles can have a vanishingly small radius), we have 
verified that the lower cutoff is sufficiently small that 
changes to its value have negligible effect on the overall 
results. 

Depending on the degree of polydispersity, the contin- 
uous parent distribution is represented by a few hundred 
species, ensuring always that the maximum sized species 
has a sufficiently small contribution to avoid spurious ef- 
fects from the upper a cutoff. For consistency, we choose 
the upper cutoff to equal that used in the simulations. 
The mean particle diameter is always taken as the refer- 
ence length scale in the calculations. 
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III. RESULTS 

Below we present our findings for the effect of the hard 
wall on the density distribution. For the purposes of 
gauging the accuracy of the DFT predictions, we have 
confined our comparison to the (representative) choice of 
the narrow top-hat distribution (c = 0.2, 5 = 0.115) and 
the Schulz distribution (z = 5, i5 = 0.407). Results for the 
wide top-hat (c = 0.7) have been obtained only via sim- 
ulation, but serve to elucidate the effect of distribution 
shape changes at a given 5. 




FIG. 2: The particle density distributions p(cr, z) obtained 
from the MC simulations. (a,b) Top-hat {5 = 0.115), (c,d) 
Schulz. 

An initial impression of the effects of confinement on 
the size distribution can be gained from fig. 2, which 
shows the form of p((7, z) obtained in the simulations for 
the narrow top hat and the Schulz distributions. Den- 
sity oscillations normal to a solid interface are a common 
feature of confined fluids, and generally increase in mag- 
nitude with increasing volume fraction. In the polydis- 
perse case, however, fig. 2 shows that the character of the 
oscillations can depend both on the particular choice of 
parent distribution and (for a given parent) on the value 
of a. 

In order to analyse the data contained in fig. 2 in 
greater detail, it is instructive to examine orthogonal 

slices p(z|cr) and p((7\z). Beginning with the former case, 
examples of p{z) for a selection of values of the particle 
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FIG. 3: (a) The spatial distribution p{z\a) for selected values 
of a for the top-hat distribution {6 = 0.115) at r?6 = 0.2. (b) 
Detail showing the comparison between MC and DFT results 
in the region close to the wall. 



diameter a are shown for the narrow top-hat distribution 
at rib = 0.2 and rib = 0.4 in figs. 3 and 4 respectively. The 
corresponding data for the Schulz distribution is shown 
in figs. 5 and 6. Included also in those figures are the pre- 
dictions of the DFT theory (cf. sec. II B). In all cases one 
observes the presence of density oscillations which start 
from the bulk density value far from the wall and grow 
with decreasing z down to the point z = a/2, where the 
curves terminate abruptly due to the excluded volume 
constraint of the hard wall. Also shown for comparison 
in each case is the form of p{z) for monodisperse hard 
spheres (with a = 1.0) at the same respective volume 
fraction and normalized so that its value at z = L/2 
matches the density profile of the a = 1.0 species in 
the polydisperse system. Additionally, we provide the 
cr- averaged particle number density {p{z\a))„ as well as 
(in cases where it deviates appreciable from the latter), 
the local volume fraction r]{z) = f / dap{z\a)a^. 

With regard to the general features of the profiles 
shown in figs, 3-6, we note that an aspect peculiar to the 
narrow top-hat parent distribution is that p{z) oscillates 
with a period a ~ 1.0 (i.e. close to a) for all particle sizes. 
This contrasts with the case of the much wider Schulz 
distribution, where the oscillations do not appear to be 
controlled by the average particle size, having a consid- 
erably larger wavelength in z. We note further that in all 
cases studied, the oscillations in p{z) induced by the wall 
are more pronounced for the monodisperse fluid than for 
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FIG. 4: (a) The spatial distribution for selected values 

of the particle size for the top-hat distribution (5 = 0.115) at 
r]h = 0.4. (b) Detail showing the comparison between MC 
and DFT results in the region close to the wall. 
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FIG. 5: (a) The spatial distribution p{z\a) for selected values 
of the particle size for the Schulz distribution at r]b = 0.2. (b) 
Detail showing the comparison between MC and DFT results 
in the region close to the wall. 



the poly disperse ones. Thus the effect of polydispersity 
appears to be to dampen the profile oscillations. This is 
presumably a result of structural disordering effects de- 
riving from the ability of small particles to occupy the 
gaps between large ones. These in turn serve to disrupt 
the excluded volume packing effects which underly the 
density oscillations commonly observed in monodisperse 
systems near a hard wall. Within this qualitative pic- 
ture, one would thus expect a greater dampening of the 
oscillations for the Schulz distribution than for the nar- 
row top-hat. on account of the disparity in their values of 
6. This indeed appears to be the case, with the top-hat 
profiles resembling more closely the monodisperse case. 

Turning next to the comparison of the simulation mea- 
surements with the results of the DFT calculations, it is 
gratifying to note that the level of agreement is generally 
excellent, particularly at the lower volume fraction. Only 
for the case rji, = 0.4 and points close to the wall do we 
observe relatively minor systematic deviations, namely a 
slightly higher local density for the MC results. 

As noted above, the profiles p{z\a) are truncated at 
z = ^ reflecting the geometrical constraint associated 
with the impenetrable wall at z = 0. The density at 
this point, pc{cr) = p{z = §|cr) is the so-called contact 
value. In fig. 7 we plot pdo') for selected a for both 
the narrow and wide top-hat, as well as the Schulz dis- 
tribution. From the figure, it is apparent that for both 
top-hat distributions, pdo') increases as a function of cr, 
i.e. the largest particles are preferentially favored over 



small ones at the wall. A similar effect occurs for the 
Schulz distribution, fig. 7(b), although in this case since 
the parent distribution is not constant in cr, one must 
form p{z = a/2\a)/p^{a) in order to expose the density 
enhancement at contact (sec fig. 7(c)). 

The contact values are related to the partial pressures 
and hence the total pressure via the sum rule [33]: 



0p = j dapdcj) (13) 

The resulting estimates of the pressure are presented 
in table I where they are compared with the predictions 
of the moment based expressions of Mansoori et al (BM- 
CSL) [34, 35] and Salacuse and Stell [8]. The former is 
based on a generalization to mixtures of the Carnahan- 
Starling equation of state, while the latter is a general- 
ization of the Percus-Ycvick (PY) theory. The simula- 
tion estimates show good agreement with the analytical 
equations of state at the lower volume fraction rjb = 0.2, 
though at the higher value there are discrepancies, with 
the BMCSL equation apparently predicting the observed 
pressure more accurately for the top-hat distribution and 
the PY equation faring better for the Schulz distribu- 
tion. The latter result is in accord with previous sim- 
ulation studies of the bulk equation of state for Schulz 
distributed hard spheres [21]. As regards the estimates 
for the pressure derived from the DFT contact values, we 
note the close agreement of these with the PY equation of 
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FIG. 6: (a) The spatial distribution p{z\a) for selected values 
of the particle size for the Schulz distribution at r]i, = 0.4. (b) 
Detail showing the comparison between MC and DFT results 
in the region close to the wall. 



state. This finding is, however, not too surprising given 
that the specific Roscnfcld functional we employ is tuned 
to recover the PY equation in the bull?; (see Sec. II B). 

The enhancement (relative to the bulk) of the larger 
particles at contact (see fig. 7(c)), is attributable to at- 
tractive depletion forces. It is interesting to note that for 
a given ijb, p{z = a/2\a)/p°{a) is apparently very similar 
for the Schulz distribution and the wider of the two top- 
hat parents. Since both these distributions have a very 
similar value of S, this suggests that the depletion inter- 
actions at the wall for a given a are primarily controlled 
by the values of rjb and S, and are much less sensitive to 
the overall distribution shape. 

The wall-fluid depletion interaction can be further 
quantified (see for example ref. [36]) in terms of an ef- 
fective potential (potential of mean force) via: 



^=exp(-/3$e//(^k)). 



(14) 



Here ^eff{z\a) represents the excess grand potential of a 
particle of size a at distance z from the wall with respect 
to its value in the bulk. Plots of ^eff{z\cr) (fig. 8) indeed 
show the expected attractive nature of the effective po- 
tential at contact, and the deepening of the well depth 
with increasing a. Note, however, that the effective po- 
tential $e// defined in eq.(14) should not be confused 
with a true two-body depletion potential. A relationship 
similar to eq.(14) has been used to quantify fluid- wall de- 



pletion effects in asymmetrical binary mixtures compris- 
ing large (colloid) particles in a solvent of small particles 
(such as polymer). There one defines a depletion poten- 
tial between large particles and the wall in the limit in 
which the concentration of the large particles [28, 36] van- 
ishes, but at fixed concentration of the smaller ones. By 
contrast, in the present case, the depletion forces are not 
mediated by a separate solvent but arise from the effects 
of the distribution of particle sizes [37] . Nevertheless sim- 
ilarities are evident with the form of the true two-body 
depletion potential observed in binary mixtures [28, 36] . 
We note further, that the potentials for the Schulz distri- 
bution (6 — 0.407) are numerically very similar to those 
for the wide top hat {S = 0.404) for distances z < 2.0. 
This suggests that our observation (see above) concerning 
the apparent insensitivity of the contact enhancements 
to the shape of /(cr), extends to the form of the effective 
potential at significant distances away from the wall. 

We turn next to an analysis of the distribution of par- 
ticle sizes as a function of the distance z from the wall, 
as expressed through the form of p{a\z). In figs. 9 and 
10, we plot p{a) at selected values of z for both the nar- 
row top-hat and Schulz parent distributions at bulk vol- 
ume fractions rji, = 0.2 and rn, = 0.4. Focusing first 
on the top-hat case, one observes that as z is increased 
(starting from near the wall), the form of p{a) alternates 
between being monotonically increasing and monotoni- 
cally decreasing. This implies that with increasing z the 
biggest and smallest particles are alternately favored, an 
effect known as local size segregation [24] (see also be- 
low). We note that at the closest allowed distance from 
the wall, large particles are always strongly favored. For 
large z, p{a) of course approaches the flat top-hat parent 
form. 

In contrast to the case of the narrow top-hat, for the 
Schulz parent the density of the largest permitted par- 
ticles (i.e. those allowed by the geometrical constraint 
at a given z) always exceeds the bulk density at that a 
(cf. fig. 10). The degree of enhancement is considerably 
stronger for r]b = 0.4 than 7^5 = 0.2. In both cases, for 
sufficiently small z, the distribution increases monotoni- 
cally. However, as z becomes larger, the enhancement in 
the density of the largest particles gradually reduces, and 
a peak (associated with the maximum of p^icr)) starts to 
appear at z ~ 1.0, signalling the crossover to bulk be- 
haviour. It is interesting to note that for % = 0.4, the 
distribution for 1.05 < 2 < CTc is very broad and exhibits 
two maxima. One of these is related to the maximum of 
the parent distribution while the other corresponds to the 
depletion interaction-induced enhancement of the density 
of the largest particles permitted at that z. Again, we 
find that for both parent forms the DFT appears to cap- 
ture semi-quantitatively the behaviour observed in the 
simulations with some relatively slight discrepancies to 
be found only at rji, = 0.4. 

In order to analyze further the local size segregation 
effects identified above, we have considered the spatial 
dependence of the local relative concentration defined as: 
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Parent 


MC 


DFT 


BMCSL 


PY 


top-hat (c — 0.2, rjt — 


0.2) 


0.8718±0.008 


0.8748 


0.8774 


0.8774 


top-hat (c = 0.2, rji, = 


0.4) 


4.955±0.074 


5.074 


4.971 


5.180 


top-hat (c = 0.7, rib = 


0.4) 


2.809±0.085 




2.86 


2.966 


Schulz (rib = 0.2) 




0.518±0.009 




0.5149 


0.5175 








0.516 


(0.5147) 


(0.5173) 


Schulz {r]b = 0.4) 




2.767±0.080 




2.668 


2.763 








2.732 


(2.593) 


(2.762) 



TABLE I: Estimates of the pressure as calculated from the contact values (eq. 13) for the MC and DFT calculations for 

the top-hat and Schulz parents at the given volume fraction r/t. Also shown for comparison are the predictions of the BMCSL 
[34] and PY [8] equations of state. The values in brackets refer to the BMCSL and PY predictions for a Schulz parent having 
a lower cutoff in a as employed in the DFT calculations (see sec. II B). 



p{z,a) 



J^da' p{z,a')' 



(15) 



This quantity measures the concentration of species of 
size (7 at a given z. Size segregation is signalled by the 
appearance of an oscillatory structure in (p{z\a). We have 
obtained the form of (/)(z, cr) for both top-hat distribu- 
tions (figs. 11,12) and the Schulz distribution (figs. 13). 
The oscillations are more pronounced at the larger rjf,. 
For the narrow top-hat, a phase difference of approx- 
imately TT radians is observed between the largest and 
smallest permitted particles and the period is close to 
the average particle size [38] . However, for the wider top 
hat, the degree of anti-correlation is somewhat less and 
the period is extended. For the Schulz parent, compar- 
ison of the curves for (/>(z|cr) is complicated by the fact 
that at large z each curve converges to a limit that is pro- 
portional to p^{(t). Nevertheless, for a given ct, clear os- 
cillations are visible, although they die out more rapidly 
than for the top hat and do not exhibit a well-defined 
anti-correlation. 

Finally in this section we turn to our results for the 
local polydispersity 5{z) defined via a generalization of 
eq. 3: 



8{z) = 



[Jda{a-afp{<y\z)\ 
J daap{a\z) 



1/2 



(16) 



Our measurements of this quantity for both top-hat 
parents and the Schulz parent are presented in fig. 14. 
One observes oscillations in 5{z), derived (cf. eq.l6) from 
oscillations in both the local mean and the standard de- 
viation of the density distribution. The oscillations are 
stronger for the wider parents than for the narrow one. 
Further insight into this behaviour can be gained from a 
perturbative analysis of 5{z) in the limit of a narrow par- 
ent. We assume that the species deviate slightly from the 
mean species, so that a = a{l + e) with e a small param- 
eter. In this limit, the density profile can be expressed 
[24] as 



= z)p° (cr) ^-l3\V{a,z)-V(&,z)] 

p{a, oo) 
X {l + e[d'{z)-c'{oo)]} 



(17) 



where p{d-, z) is the density profile for the mean species, 
while c' corresponds to the reversible work associated 
with slightly changing the size of one particle, a distance 
z from the wall, within an otherwise monodisperse sys- 
tems at the corresponding volume fraction (see [24] for 
further details). Introducing this expansion into eq.(16) 
yields 



6{z) 



I p{a, oo)e2 
p{a,z)no 





e3 — 


X |l + 


2e2 



[c'{z) - c'(oo) 



(18) 



at distances where the wall potential is negligible. 

Eq. 18 indicates that for a narrow parent the local 
deviations in the width of the distribution depend on 
second and higher moments of the parent. As a re- 
sult, in this limit the local variations of the mean density 
make a relevant contribution to 6. This is confirmed by 
fig. 14(a) which shows that the oscillations are indeed 
anticorrelated with those of the mean species. In fact, 
the amplitude of the term c'{z) — c'(oo) is around 0.004 
for a flat parent with 5 = 0.115 and around 0.054 for 
S = 0.404. Since the deviations of c' are never larger 
than unity for the parameters considered (see e.g. Fig. 
2 in ref.[24]), this explains the relative weakness of the 
oscillations in fig. 14(a). It also helps to explain the be- 
havior of fig. 14(b) (although in this case one is far from 

2 

the narrow limit): for a Schulz distribution = 26^ and 

hence the correction in eq.(18) vanishes. Accordingly, in 
a Schulz distribution the shape corrections to 6{z) will 
become significant only at large S. In fact, comparing 
S{z) for S = 0.404 between the Schulz and top hat par- 
ents (fig. 14), one can see that for a Schulz the shape of 
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FIG. 7: Contact values Pc{cr) = p{(7,z = cr/2). (a) Top-hat 
distributions {S = 0.115 and S = 0.404) at volume fractions 
rib = 0.2,0.4. (b) Schulz distribution {5 = 0.407). (c) The 

contact values relative to the bulk density. In all cases, lines 
are merely guides to the eye. Error bars are smaller than the 
symbol sizes. 



the local degree of polydispersity follows more closely the 
density profile of the mean species. 

Although wc have used a definition of the local poly- 
dispersity that generalized the conventional form (eq.3) 
used to describe bulk parent forms, it is interesting to 
note that, at least for inhomogeneous systems, a defi- 
nition of d{z) based on the local concentration, (/)(cr, r), 
rather than on the local number density, serves to better 
illustrate the role of polydispersity effects. Accordingly 
we form a narrow limit expansion of 



S{z) 



[J da{a — (t)^(I){cf, z)] 
J daa<p{a, z) 



1/2 



(19) 



FIG. 8: The effective wall-particle potential at volume frac- 
tion rib = 0.4. (a) Narrow top-hat (<5 = 0.115); (b) Schulz, 
(5 = 0.407) and (c) Wide top-hat {5 = 0.404) distribution. 



which yields, 





" £3 — 




2 




2? 



[c'(z)-c'(oo)] . (20) 



One sees that this last expression does not depend on the 
local density of the mean species; accordingly we would 
generically expect a much weaker variation of 6(z) for a 
Schulz, and for other asymmetric parents. This alterna- 
tive definition of the local polydispersity thus highlights 
the fact that in the narrow limit the local segregation of 
species tends to preserve the width of the parent imposed 
in the bulk. 



IV. DISCUSSION AND CONCLUSIONS 

In summary we have employed specialized MC simula- 
tions techniques and DFT to study the effects of a hard 
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FIG. 9: The size distribution at selected distances z from 
the wall for the narrow top-hat distribution (5 = 0.115). (a) 
rit = 0.2, (b) rjt = 0.4. Symbols: Monte Carlo results. Solid 
lines: DFT results. 



wall on the properties of polydisperse hard spheres. Bulk 
phase (parent) distributions of the top-hat and Schulz 
forms were considered. In the former case, degree of poly- 
dispersity S = 0.115 and 6 = 0.404 were studied, while 
in the latter case, S = 0.407 was used. In all cases the 
density distribution p{a, z) was obtained and analyzed 
for bulk volume fractions r/(, = 0.2 and % = 0.4. 

The original motivation for this study was a desire to 
gauge the accuracy of the DFT predictions of ref. [24] via 
a like-for-like comparison with simulation. The results 
presented in sec. Ill demonstrate that the agreement is 
extremely good. Indeed, the DFT calculations provide a 
quantitatively accurate description of the system prop- 
erties at ryb = 0.2 for both parent forms studied. Even 
for the higher volume fraction r/ft = 0.4, the agreement is 
semi-quantitative. Clearly this finding bodes well for the 
future utility of DFT in investigating other polydisper- 
sity related issues in inhomogeneous fluid. 

Beyond this, the present study extends that of ref. [24] 
by considering the nature of the attractive depletion in- 
teractions in the vicinity of the hard wall. Analysis of 
contact values showed a clear preference for the largest 
particles to occupy positions in which their surfaces touch 
the wall. Indeed for the Schulz distribution, the contact 
values of particles at the cutoff exceeded that in the bulk 
by over two orders of magnitude. Interestingly, at a given 
77fc, the degree of contact enhancement as a function of 
cr, seems quite insensitive to the distribution shape (cf. 
fig. 7(c)) for a given 5. This finding was found to ex- 
tend to the form of the effective potentials (fig. 8) near 
the wall. Clearly however our comparison is far from 
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FIG. 10: The size distribution p{a\z) at selected distances z 
from the wall for the Schulz distribution {5 = 0.407). (a) 
77b = 0.2, (b) rjb = 0.4. Symbols: Monte Carlo results. Solid 
lines: DFT results. 

exhaustive in this regard and further studies would be 
required to determine whether this observation applies 
more generally. 

The interplay of the cr-dependent attractive depletion 

forces and the geometrical constraint imposed by the wall 
was found to radically alter the local size distribution. 
Indeed in the case of the Schulz parent at % = 0.4, a 
minimum in p(cr\z) was observed to develop within a cer- 
tain range of z. Another interesting feature, was the 
appearance of oscillatory structure in the local degree of 
polydispersity ^{z). This was accompanied, very close to 
the wall, by considerable reduction in the local degree of 
polydispersity with respect to the bulk. It is tempting 
to speculate that this latter finding may have some bear- 
ing on recent experimental observations of wall induced 
freezing in a polydisperse system of colloidal hard spheres 
that apparently forms a glass in the bulk [39]. Polydis- 
persity is known to hinder freezing in colloidal systems 
[6], and consequently any narrowing of picf) due to the 
presence of a wall may serve to promote the formation of 
an ordered crystal phase. 

With regard to the size segregation effects originally 
reported in ref. [24], we find that the anti-correlation in 
the oscillations of the local concentration between the 
largest and smallest permitted particles is significantly 
more pronounced for small b than for large 5 at a given 
ryt. Furthermore as for small (5, the period of the oscilla- 
tions in z appears to be close to that of the mean parti- 
cle diameter. Both these findings are in accord with the 
predictions of perturbative calculations [24]. For large 
(5, by contrast, the present work shows that distribution 



11 





FIG. 11: The local relative concentration <j){z\a) for a nar- FIG- 13: The local relative concentration (/)(2;|o-) for the Schulz 

row top-hat parent distribution [S = 0.115). Data are shown parent {S = 0.407). Solid lines are the DFT predictions, (a) 

for a selection of particle sizes a. Solid lines are the DFT rib = 0.2; (b) r)b = 0.4. 
predictions, (a) r)b = 0.2; (b) ijb = 0.4. 
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shape effects are important in determining the character 
of size segregation effects. Thus for the Schulz parent at 
S PS 0.4, the degree of anti-correlation is much less than 
for the corresponding top-hat parent. Additionally, in 
the latter case the period of the oscillations appear to 
be more consistent with the average particle size than 
for the Schulz distribution. This suggests that it is the 
presence in the Schulz parent of relatively large higher or- 
der moments which leads to this increased disruption of 
size-segregation effects. Analogous observations apply to 
the density profiles p{z\(j) themselves. In all cases, (but 
especially for the Schulz parent), the presence of poly- 
dispersity was found to dampen the density oscillations 
significantly compared to the monodisperse profile p{z). 

Turning finally to the outlook for further related work. 



FIG. 14: The spatial dependence of the local degree of poly- 
dispersity S{z) at volume fraction rit = 0.4. (a) Top-hat, 
S = 0.115; (b) Schulz S = 0.407; (c) top-hat, S = 0.404. 
Also shown for the respective parents ((a') etc) are the mean 
(circles) and the standard deviation (squares) of the local size 
distribution p{a\z). Lines show the profile of the mean species 
p{z\a = 1), relative to its bulk value. 



an interesting open question is whether size segregation 
effects analogous to those considered here for an inhomo- 

geneous system might also occur in bulk polydisperse flu- 
ids i.e. whether pairs of hard spheres of prescribed sizes a 
and a' are preferentially found with certain separations. 
This could be investigated by appeal to the measured 
form of the pair distribution function g{r, a, a'). We hope 
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to report on this matter in a future publication. 
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